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Abstract: 

A previous study analyzed errors in the numerical calculation of actual crop evapotranspiration (ET a ) under soil water stress. 
Assuming no irrigation or precipitation, it constructed equations for ET a over limited soil-water ranges in a root zone drying out 
due to evapotranspiration. It then used a single crop-soil composite to provide recommendations about the appropriate usage of 
numerical methods under different values of the time step and the maximum crop evapotranspiration (ET C ). This comment 
reformulates those ET a equations for applicability over the full range of soil water values, revealing a dependence of the relative 
error in numerical ET a on the initial soil water that was not seen in the previous study. It is shown that the recommendations 
based on a single crop-soil composite can be invalid for other crop-soil composites. Finally, a consideration of the numerical 
error in the time-cumulative value of ET a is discussed besides the existing consideration of that error over individual time steps as 
done in the previous study. This cumulative ET a is more relevant to the final crop yield. Published 2014. This article is a U.S. 
Government work and is in the public domain in the USA. 
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THE SHANG (2012) STUDY: ASSUMPTIONS, 
EQUATIONS AND RECOMMENDATIONS 

Shang (2012) is hereforth denoted in this comment article 
as S12, whereas this article itself is denoted as YC for 
Yatheendradas (et al.) comment. Besides the original 
assumptions for calculating the actual crop evapotranspi- 
ration (ET a ) by Allen et al. (1998), S12 made important 
additional assumptions that are accordingly followed here 
as well: 

1 . Meteorological/atmospheric conditions were constant 
over multiple time steps (see Figure 3 by S12). This 
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means that the standard reference crop evapotranspiration 
(ET 0 ) from either grass or alfalfa was assumed constant. 

2. Important source terms were neglected in the 
analytical and numerical formulations, where their 
inclusion may change the results and recommenda- 
tions. One such term is the water input (irrigation or 
precipitation). This essentially slowly dried out the 
crop root zone because of evapotranspiration, 
starting from whichever considered instant in the 
crop phenology. 

3. Another source term or key parameterization not 
considered but existing in reality and in practice is 
the time variation in the crop coefficient curve (Allen 
et al., 1998). Typically, this curve assumes a low 
constant initial value phase, followed by a linear 
increase during a crop development phase, then a 
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maximum constant value over a mid-season phase, and 
finally a linear decrease till harvest in a late-season 
phase. The instantaneous maximum crop evapotrans- 
piration (ET C ) is obtained by multiplying ET 0 by this 
crop coefficient. So for a constant ET 0 from assump- 
tion 1 aforementioned, a constant ET C over multiple 
time steps implies a corresponding enforced constancy 
in the crop coefficient. Figures 1-3 by S 12 are obtained 
using this enforcement. 

4. S12 neglected soil water flux through the bottom of 
root zone by assuming a ‘deep groundwater level and 
in periods without heavy rain and irrigation’. 

Following Allen et al. (1998), Equation (2) by S12 
gave ET a using ET C and a water stress coefficient (K s ): 

ET a - K S ET C (1) 


Wj = pW p + (1 — p)W f , 0 < p < 1 (3) 

where p is the depletable fraction of root zone water 
before water stress, and Equation (5) by S12 showed p 
depending both on the crop type through a p std term and 
on the ET C (mm/day): 

p - Pstd + 0.04(5 - ET C ) (4) 

Finally, assuming that soil water stayed in the 
W p < W < Wj range in Equation (2) earlier for a time 
step. Equation (8) by S12 gave the following analytical 
solution (ET b ) of ET a : 


ET b = (W 0 - W, 



r 

ET C At 

1 

Wp)< 

r e 

Wj-WpJ 



>,W p < W< Wj (5) 


Equation (3) by S12 gave this piecewise linear K s that 
increases with an increase in the soil water storage (W) 
from the wilting point (W p ) to the critical soil water 
storage for soil water stress (Wj): 

f 0, W<W P 

K s = < (W - W p )/(Wj - W p ), W p < W < Wj (2) 
( 1, W>Wj 


where W Q is the initial soil water storage, At is the time 
step, and ET C is the ET C assumed constant (or 
‘representative’) over the time step. Note that although 
ET C is a rate, its analytical solution ET b was formulated 
as a depth (mm). 

For the same W p < W < Wj range, the Equations (9), (15) 
and (17) by S 12 gave the following numerical solutions of 
the ET a depth (mm), denoted by ET ee , ET me and ET H 3 for 
the explicit Euler, the modified Euler and the Heun’s third- 
order methods, respectively: 


Flere, Wj was given by Equation (4) by S12 as lying 
between W p and the soil water at field capacity (Wf): 


ET ee 


w 0 - Wp 

Wj - Wp 


ET C At, 


W p < W < Wj 


( 6 ) 


-nit— S12, Wo=200 mm 
-his— S12, Wo=220 mm 
—sis— SI 2, Wo=240 mm 
— *— SI 2, Wo=260 mm 


(a) Pstd=°-55 



-□- YC, Wo=200 mm 
-E- YC, Wo=220 mm 
-E- YC, Wo=240 mm 
- E- YC, Wo=260 mm 


(b) p std =0.8 



Figure 1. Analytical ET a (i.e. ETb) for different values of ET C , W Q and p st d. In the legend, S12 and YC, respectively denote usage of equations from 

Shang (2012) and this comment article 
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SI 2: ET =5 mm/d i— i YC: ET =5 mm/d 
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SI 2: ET=7 mm/d YC: ET=7 mm/d 


SI 2: ET =9 mm/d 


(a) Explicit Euler method: 
p =0.55, W =200 mm 


YC: ET =9 mm/d 

C 

(b) Explicit Euler method: 
p =0.55, W =260 mm 

'std o 



Figure 2. Relative error in ET a for some numerical methods, crop-soil composites, initial conditions and ET C values. In the legend, S12 and YC, 
respectively, denote the usage of equations from Shang (2012) and this comment article. Horizontal dashed line is the relative error threshold at 5% 

(or —5% as relevant) 
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( 8 ) 


S12 further used analytical Equation (5) in the previous 
texts as a benchmark to calculate the relative errors in the 
numerical ET a from these Equations (6)-(8). According- 
ly, S12 preferred less complex numerical methods for 


computational efficiency in making the following recom- 
mendations on their selection for calculating a satisfac- 
torily accurate numerical ET a : 

A. If the time step is 1 day (d) and ET C is any value, then 
use explicit Euler, else, 

B . if the time step is up to 2 days and ET C is below 5 mm/day , 
then use explicit Euler, else, 

C. if the time step is up to 1 week and ET C is any value, 
then use midpoint/modified Euler, else, 

D. if the time step is up to 10 days and ET C is below 
5 mm/day, then use midpoint/modified Euler, else, 

E. if the time step is up to 1 week and ET C is any value, 
then use Heun’s 3 rd -order method. 
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. W o =200 mm, 1-day _ 0 _ W^-260 mm, 1-day 

W o =200 mm, 2-day W o =260 mm, 2-day 

W =200 mm, 1 -week i— i W =260 mm, 1 -week 

o □ o 

W o =200 mm, 10-day W o =260 mm, 10-day 

W o =200 mm, 15-day _ 0 _ W o =260 mm, 15-day 


(a) Explicit Euler, ET C = 5 mm/d 


(b) Explicit Euler, ET C = 9 mm/d 



20 40 e 

(c) Explicit Euler, ET C = 1->9 mm/d 


20 40 60 

(d) Explicit Euler, ET C = 9->1 mm/d 



Figure 3. Relative error in cumulated ET a of S12’s crop-soil composite using the reformulated explicit Euler method for different ETc trends, initial 
conditions (Wo) and time step. Horizontal dashed lines are the relative error thresholds at ±5% 


These recommendations were obviously made with 
practical applicability in mind. Hence, it is worth noting 
that following them for computational efficiency while 
attempting a realistic case of a non-constant ET C over 
multiple time steps (i.e. by relaxing assumption 3 earlier) 
can mean continuously changing the numerical method 
during the simulation. For example, for a fixed time step 
of 2 days, a change in ET C from 4 mm/day at a time step 
to 6 mm/day at the next would correspond to a change in 
the numerical method from the explicit Euler to the 
midpoint/modified Euler. 

REFORMULATING THE S12 EQUATIONS: EFFECT 
ON RESULTS AND RECOMMENDATIONS 

Figure 1 by S12 showed the analytical ET b after a 10-day 
time step over a range of values for the ET C (1 to 9 mm/ 
day) and initial soil water (W Q from 200 to 260 mm). 
These ranges and the properties for crop (p std = 0.55) and 
soil (W p = 158 mm, Wj = 396mm) are characteristic of 
their study site in North China. However, some 
combinations of ET C and W 0 actually move W 0 out of 
the assumed W p < W < Wj range and into the W > Wj 
range. For example, an ET C of 1 mm/day gives a Wj of 
around 227 mm using Equations (3) and (4) earlier. 


compared with which a W 0 of 260 mm is greater and so 
outside the W p < W < Wj range. Such combinations of 
ET C and W D are omitted in Figure 1 by S12 and also in 
Figure la here wherein some asterisks derived using 
Equation (5) earlier are conspicuously absent at some 
square marker locations (these square markers are based 
on a reformulated ET b to be presented further in the 
succeeding texts). These combinations of ET C and W 0 
having W Q > Wj tend to occur on the higher W Q curves at 
lower ET C values. Note that at higher ET C values on the 
same W 0 curves, W p < W 0 < Wj happens because of 
Equations (3) and (4) earlier that result in a Wj that 
increases with increasing ET C . 

For crop or soil properties different from the specific 
values considered by S12, such missing ET b values due to 
the W p < W < Wj range constraint can occur for a 
significantly large number of combinations of ET C and 
W D . For example, see asterisks absent at more than half of 
the square marker locations in Figure lb here, wherein the 
crop p std value is now a higher-end value of 0.8 (Allen 
et al., 1998). 

Analytical and numerical formulations of ET a specific 
to limited soil water ranges or not accounting for 
transition between such ranges would be limited in their 
applicability. For example, a formulation valid for soil 
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water in the W > Wj range can be applied only if both W Q 
and the soil water at the end of a time step are in that 
range. Avoiding this limited applicability requires 
inclusion in the formulations of the transition between 
ranges by soil water during a time step, for example, 
from the W > Wj range to the W p < W < Wj range. Also, 
the ET b calculated using Equation (5) earlier for a W Q in 
the W > Wj range can be overbiased near ET C -W 0 
combinations where such transition occurs. As a 


At, = (W 0 - Wj)/ ET C , Wo > Wj (9) 

Note that soil water hitting the lower threshold W p 
during a time step is handled by enforcing this W p 
constraint using ‘min' functions. 

Applying this time step and ET a splitting and the W p 
constraint gives the following generally applicable analytical 
solution (ET h c in mm) of ET a instead of the Equation (5) 
earlier by S12 that had a W p < W < Wj applicability only: 


r o. 


Wo<W p 


ET b c = < 


mm 


mm 


ET c At 


(W 0 - W p ) / 1 - eL w j - w rJ Uw 0 - W p ) 


ET C At, + (Wj-WpW 1 - eL 


ET c (At — At,) 


Wj-W p 


,W P < W 0 < Wj 


(W 0 - W p 


ET C At, 


(W 0 >Wj)&([W 0 -ET c At] <Wj) 
(W 0 >Wj)&([W 0 -ET c At] >Wj) 


( 10 ) 


benchmark to calculate the relative error in numerical ET a , 
accurate ET b is important. Hence, the analytical and 
numerical S 12 equations are reformulated in the succeeding 
texts to avoid the aforementioned applicability constraints. 

For a time step over which the soil water is above Wj 
throughout, the ET a is constant at ET C . For a time step that 
starts and ends with soil water values above and below Wj, 
respectively, the simpler ET a equations due to the S12 
assumptions enable splitting the ET a into two terms. The 
initial term’s corresponding portion of that time step when 
soil water is above Wj is: 


These ET b c values are the ones mentioned further 
earlier as being plotted in Figure 1 here using square 
markers. When W p < W < Wj over a time step, they 
coincide with asterisk markers derived using Equation (5) 
earlier by SI 2. 

Similarly, the following are the numerical solutions of 
ET a in millimetre (denoted by ET ee c , ET me c and ET ro , 
c when using the explicit Euler, the modified Euler 
method and the Heun’s third-order methods, respectively), 
instead of the Equations (6-8) earlier by S12 that had a 
W p < W < Wj applicability only: 


ET 


ee,C — 


o, 
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'W c - W p 
Wj - Wp 


ET C At, (W 0 



W 0 <W p 
Wp < W 0 < Wj 


min [ET C At, (W„ - W p )] , (W 0 >Wj)&([W 0 - ET C At] < Wj) 
k ET C At, (W 0 >Wj)&( [Wo - ET C At] >Wj) 


( 11 ) 


Published 2014. This article is a U.S. Government work and is in the public domain in the USA. 


Hydrol. Process. 28 , 3833-3840 (2014) 


3838 


S. YATHF.F, NOR ADAS ET AL. 


f 0, 


ET me ,c = < 


Wo Wp ET r At( 1 - 1 ) , (W 0 - W p ) 


Wj -W p 


2 Wj — W p 


w 0 <w p 

,W p < W 0 < Wj 


ET r 
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Note that in Equation (11) earlier, the apparently missing 
At, term for the [(W 0 >Wj)&( [W 0 — ET c At] < Wj)] case 
is because the original expression [ET c At] + ET c (At — At, )] 
reduces to ET c At. 

When formulating the relative error in the numerical 
ET a , S12 noticed that the common (W 0 -W p ) term in the 
analytical Equation (5) earlier and the numerical Equa- 
tions (6-8) earlier cancelled out. S12 duly asserted that 
‘relative errors of ET a calculated with numerical methods 
in one time step arc independent of the initial soil water 
storage in the range of soil water stress’, and any mention 
of the initial soil water (W D ) is correspondingly absent in 
the S12 recommendations. However, this supposed 
independence from W Q is erroneous as seen from the 
analytical Equation (10) and numerical Equations (11-13) 
earlier, where the (W 0 -W p ) term need not cancel out 
while formulating this relative error This is also 
illustrated by the example contrast of this reformulation- 
derived relative error (denoted by square markers) 
between Figure 2a (W Q = 200 mm) and 2b (W Q = 260 mm) 
here that use the explicit Euler method, for the lower ET C 
curves. Also, asterisk markers in Figure 2a here plotted 
using the S12 equations are same as the markers in 


Figure 2a by S12; it is evident that the latter plot should 
actually have been impossible with some W Q values, 
because an example W Q = 260 mm should reveal missing 
markers for the lower ET C curves as shown by the missing 
asterisk markers in Figure 2b here. This was missed out 
on by the S 12 study because its relative error formulation 
that was independent of W D did not reveal any missing 
markers 

In Figure 2c-2d here, the original field capacity (W f ) and 
wilting point (W p ) values for a soil type closer to clay as 
considered by S 12 are halved towards a soil type closer to 
loam (Walker, 1989), with W Q also correspondingly halved; 
this effectively decreases the soil water storage capacity and 
hence the water availability for evapotranspiration. The 
relative error for the explicit Euler method is now seen to 
decrease after some time step (reflecting the numerical water 
content value reaching W p ) and decreases earlier with an 
increase in p std (0.55 for Figure 2c vs 0.8 for Figure 2d). This 
contradicts the trend from the equations and Figure 2 of the 
S12 study where the relative error seemed to monotonically 
increase with an increasing time step. 

The high p std value of 0.8 in Figure 2d here collapses 
the relative error value curves for ET C of 5-9 mm/day into 
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a single one (not visually discernible) that is distinct from 
the ET C curves of 1-3 mm/day. This also happens for the 
example plots for the modihed Euler (Figure 2e) and the 
Heun’s third-order methods (Figure 2f) having this same 
Pstd- Note that the usage of the original equations by SI 2, 
where occurrence of soil moisture values outside the valid 
range were ignored, would have collapsed all the ET C 
curves of 1-9 mm/day into this single curve. Hence, the 
reformulated equations accurately and better resolve the 
behaviour in relative error for such high p std values. 

A comparison between Figure 2b-2d here shows the 
relative error curves starting out steeper for lower soil 
water storage capacities (i.e. lower W f and W p ) and 
higher p std values. This aspect was not highlighted by 
S12; although the properties of the specihc crop-soil 
composite used were specified, the absence of a caveat 
about the validity of those results regarding relative error 
for that specific composite clearly gave an erroneous 
impression that those results and consequent recommen- 
dations were valid for all crop-soil composites. 

Using the same relative error threshold of 5% by S12 
and a time step of 1 day, Figure 2c-2d here show that the 
explicit Euler method can only be used for these crop-soil 
composites if ET C is some respective specific value less 
than 5 mm/day. This contrasts with the recommendation 
‘A’ earlier by S12 that this method can be used for any 
ET C for a 1-day time step. Similarly, for a 2-day time step, 
the explicit Euler method is supposedly invalid if ET C is 
above 5 mm/day as per the recommendation B' above by 
S12, but Figure 2c-2d here show this method becoming 
invalid above some lower respective specific value of ET C 
between 3 and 5 mm/day and dependent on the crop-soil 
composite and its initial condition. 

Similarly, for the modified Euler method, the respective 
recommendations ‘C’ and ‘D’ earlier by S12 state that the 
method can be respectively used for any ET C in case of a 
1 -week time step and for an ET C of 5 mm/day or less in 
case of a 10-day time step. However, for the specific crop- 
soil composite and initial condition in Figure 2e here, this 
method can be used only when ET C is below some 
respective specific value between 1 and 3 mm/day for 
both time steps. Finally, the recommendation ‘E’ earlier 
by S12 states that the Heun’s 3 rd -order method can be 
used for up to a 15-day time step for any ET C . However, 
this recommendation becomes invalid when ET C is above 
3 mm/day for the specific crop-soil composite and initial 
condition in Figure 2f here (note that its x-axis extends 
out to 15 days). 

CONSIDERING THE RELATIVE ERROR IN 
CUMULATIVE ET A 

The cumulated ET a at the end of the crop season (or harvest) 
is more relevant to the crop yield than the ET a over any 


single time step and is a common metric used in drought 
early warning systems (e.g., Senay and Verdin, 2002). In 
providing recommendations, S12 considered only single 
time-step simulations, with the largest time step being 
15 days. For considering multiple time-step simulations, a 
simulation period of 60 days (approximately 2 months) is 
used here. 

Crop seasons are typically quite longer than 60 days, with 
some individual crop growth phases possibly spanning up to 
60 days or even more (Allen etal., 1998). The assumption 3 
earlier of a constant ET C even over multiple time steps as 
shown in Figure 3 by S 12 is now relaxed here to allow a 
different ET C for each time step in a simulation. This allows 
a piecewise linear increase/decrease of ET C over a 
simulation and correspondingly a characterization of 
linearly increasing/decreasing portions of the crop coeffi- 
cient curve. Any desired temporal combination of these 
portions with constant ET C portions enables an assessment 
of the relative error in the cumulated ET a over the relevant 
crop phenology simulation. 

For the specific crop-soil composite used by S12, the 
60-day evolution of the relative error in cumulated ET a 
using the reformulated equations for the explicit Euler 
method is shown here for the following: a constant ET C at 
5 mm/day appearing in the S12 recommendations 
(Figure 3a), a constant ET C at constant upper limit of 
9 mm/day considered by S12 (Figure 3b), an ET C linearly 
increasing from 1 mm/day to 9 mm/day during the simula- 
tion (Figure 3c), and an ET C linearly decreasing from 9 mm/ 
day to 1 mm/day (Figure 3d). Figures 3a-3d show that each 
relative error is actually near or below 5% at 60 days for a 
time step of up to 10 days with a projected decreasing trend 
beyond 60 days, meaning that the relative error for their 
piecewise linear combination over time should also be 
within 5% at 60 days and beyond. Also note the relative 
eixor at 60 days for an ET C of 9 mm/day is actually lower 
than that for an ET C of 5 mm/day, due to a faster decrease in 
its relative error over time (Figures 3a-3b). 

SUMMARY 

The S12 assumptions included (i) ignoring water input 
resulting in a root zone that eventually dries out as with 
rain-fed crops in water-scarce regions, (ii) constant ET C 
over multiple time steps and (iii) deep groundwater 
implying negligible soil water flux through the bottom of 
root zone. This latter flux is significant in many regions 
(e.g. Schmid and Hanson, 2009; Schmid et al., 2009). 
The assumption of no water input is specifically 
restrictive against a practical application of its equations 
for estimating ET a . This essentially means that this model 
is inapplicable when water input becomes available. The 
S12 analytical and numerical formulations of ET a over 
limited soil water ranges without accounting for transition 
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between the ranges restrict their applicability and accuracy. 
Also, the erroneous impression of a general validity of the 
S12 recommendations could have been avoided through 
exercising more rigour in either checking results for other 
crop-soil composites or specifying relevant caveats. 

Although conforming to the S12 assumptions, a 
reformulation of its equations here towards applicability 
over the entire soil water range reveals a dependence of the 
relative error in numerical ET a on the initial soil water. This 
relative error is also shown to be greatly dependent on the 
crop-soil composite. The consequently more complicated 
actual nature of this relative error now shows that crop and 
soil properties and the initial soil water also have to be 
considered in making recommendations regarding an 
appropriate numerical method. 

Finally, any such recommendation can also consider 
the simulation setup and goal. For example, instead of 
considering the error in numerical ET a at every time step, 
one can alternatively consider the error in cumulative ET a 
for calculating crop yield. 
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